% A5_tinterp_mm5.m 9/27/2011 Parker MacCready
%
% this takes the regridded mm5 files and...
% fills in missing days by interpolating in time

clear

% *** USER DEFINED VALUES **********************
dth = 120; % time step for forcing fields (hours) [120, 12, 3, 1]
year_vec = 2001:2009; % span of years to regrid (could be just one)
%year_vec = 2001:2005;
%year_vec = 2006:2009;
% *** END USER DEFINED VALUES ******************

switch dth
    case 120
        indir0 = 'regridded_five_day/';
    case 12
        indir0 = 'regridded_half_day/';
    case 3
        indir0 = 'regridded_three_hour/';
    case 1
        indir0 = 'regridded_hourly/';
    otherwise
        disp('that dth is not coded for')
        break
end

% create the output directory
outdir0 = strrep(indir0,'regridded','tinterp');
if ~(exist(outdir0,'dir')==7); mkdir(outdir0); end;

disp(' '); disp(['*** Interpolating from ',indir0,' ***']); disp(' ');

% LOOPS

% choose a variable
varlist = {'psfc';'t2';'qair';'rain';'swdown';'lwdown';'u10r';'v10r'};
for vv = 1:length(varlist);
    var = varlist{vv};
    
    disp([' - working on variable = ',var]);
    tic
    
    td = []; % initialize datenum vector
    mo = []; % initialize a month vector
    yr = []; % initialize a year vector
    ind = []; % index within the month
    data = []; % initialize a data vector
    
    % choose a year
    for year = year_vec
        
        indir = [indir0,num2str(year),'/'];
        
        outdir = [outdir0,num2str(year),'/'];
        if ~(exist(outdir,'dir')==7); mkdir(outdir); end;
        
        disp(['  -- year = ',num2str(year)]);
        
        % choose a month
        for month = 1:12
            
            load([indir,var,'_',num2str(month),'.mat']);
            
            % need to make copies so that we don't overwrite the
            % original regridded files (which take a long time to make)
            save([outdir,var,'_',num2str(month),'.mat'], ...
                'VAR','LON','LAT','TD','flag_interpolated','-v7.3');
            
            % and gather some information used in the time interpolation
            nt = length(TD);
            td = [td; TD];
            mo = [mo; month*ones(nt,1)];
            yr = [yr; year*ones(nt,1)];
            ind = [ind; [1:nt]'];
            data = [data; squeeze(VAR(:,10,10))];
            
            [NR,NC] = size(LON);
            
        end % end of month loop
        
    end % end of year loop
    
    % find the gaps
    isbad = isnan(data);
    
    % interpolate across gaps
    starts = find(diff(isnan(data))==1);
    ends = find(diff(isnan(data))==-1) + 1;
    % and find the delta_t (hours)
    delta_t = round(24*(td(2) - td(1)));
    
    for ii = 1:length(starts)
                
        i0 = starts(ii); i2 = ends(ii);
        
        ninterp = i2 - i0 - 1;
        
        for jj = 1:ninterp;
            
            i1 = i0 + jj; % index of the point we are interpolating to
            
            % need to choose start and end points that are at the same hour
            % of the day as the point we are interpolating to.
            
            % first find the hours of the point we are working on (hi1) and
            % also the hours of the first good point before (hi0) and after
            % (hi2)
            hi0 = datevec(td(i0)); hi0 = hi0(4);
            hi1 = datevec(td(i1)); hi1 = hi1(4);
            hi2 = datevec(td(i2)); hi2 = hi2(4);
            
            % now make the index for the first point before that has the
            % same hour as hi1
            ii0 = i0 - abs((hi1 - hi0)/delta_t);
            % and then keep marching back by days untl we find the index of
            % a good point.  The "ceil" call is meant to handle the case
            % where delta_t is longer than a single day (e.g. 5 days)
            while isbad(ii0)
                ii0 = ii0 - ceil(24/delta_t);
            end
            % same procedure for the end point of the interpolation
            ii2 = i2 + abs((hi2 - hi1)/delta_t);
            while isbad(ii2)
                ii2 = ii2 + ceil(24/delta_t);
            end
            
            new_ninterp = ii2 - ii0 - 1;
            nninterp = new_ninterp + 1;
            
            % load the starting data slice(s)
            A = load([outdir0,num2str(yr(ii0)),'/', ...
                var,'_',num2str(mo(ii0)),'.mat']);
            AA = squeeze(A.VAR(ind(ii0),:,:));
            % load the ending data slice
            C = load([outdir0,num2str(yr(ii2)),'/', ...
                var,'_',num2str(mo(ii2)),'.mat']);
            CC = squeeze(C.VAR(ind(ii2),:,:));
            
            % for times between:
            % (i) interpolate to make the new slice
            % (ii) figure out where it goes and put it there
            B = load([outdir0,num2str(yr(i1)),'/', ...
                var,'_',num2str(mo(i1)),'.mat']);
            B.VAR(ind(i1),:,:) = ((i1-ii0)*(CC-AA)/nninterp) + AA;
            B.flag_interpolated(ind(i1)) = 1;
            % save the results
            VAR = B.VAR; LON = B.LON; LAT = B.LAT; TD = B.TD;
            flag_interpolated = B.flag_interpolated;
            save([outdir0,num2str(yr(i1)),'/', ...
                var,'_',num2str(mo(i1)),'.mat'], ...
                'VAR','LON','LAT','TD','flag_interpolated','-v7.3');
            
        end % end of the jj loop
        
    end % end of the starts-ends loop
    
    dt_toc = toc;
    disp([' - interpolation of ',var,' took ', ...
        num2str(round(dt_toc)),' sec'])
    disp('***********************************')
    
end % end of variable loop






